{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example: Static inverse free-boundary equilibrium calculations (in ITER)\n", "\n", "---\n", "\n", "Here we will generate an equilibrium (find coil currents with the inverse solver) in an ITER-like tokamak. \n", "\n", "The machine description comes from files located [here](https://github.com/ProjectTorreyPines/FUSE.jl).\n", "\n", "The equilbirium\\profile parameters are **completely made up** - please experiment on your own and change them to more realistic values as you please!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Import packages" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Create the machine object" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/ITER/ITER_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/ITER/ITER_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/ITER/ITER_limiter.pickle\",\n", " wall_path=f\"../machine_configs/ITER/ITER_wall.pickle\",\n", ")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# plot the machine\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "plt.tight_layout()\n", "\n", "tokamak.plot(axis=ax1, show=False)\n", "ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", "\n", "ax1.grid(alpha=0.5)\n", "ax1.set_aspect('equal')\n", "ax1.set_xlim(1.1, 12.4)\n", "ax1.set_ylim(-8.3, 8.3)\n", "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", "ax1.set_ylabel(r'Height, $Z$ [m]')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate an equilibrium" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import equilibrium_update\n", "\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", " Rmin=3.2, Rmax=8.8, # radial range\n", " Zmin=-5, Zmax=5, # vertical range\n", " nx=129, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Instantiate a profile object" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# initialise the profiles\n", "from freegsnke.jtor_update import ConstrainBetapIp\n", "profiles = ConstrainBetapIp(\n", " eq=eq, # equilibrium object\n", " betap=0.15, # poloidal beta\n", " Ip=11e6, # plasma current\n", " fvac=0.5, # fvac = rB_{tor}\n", " alpha_m=2.0, # profile function parameter\n", " alpha_n=1.0 # profile function parameter\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load the static nonlinear solver" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Constraints" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Rx = 5.02 # X-point radius\n", "Zx = -3.23 # X-point height\n", "\n", "Ro = 6.34 # X-point radius\n", "Zo = 0.66 # X-point height\n", "\n", "# set desired null_points locations\n", "# this can include X-point and O-point locations\n", "null_points = [[Rx, Ro], [Zx, Zo]]\n", "\n", "# set desired isoflux constraints with format \n", "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", "# with each isoflux_i = [R_coords, Z_coords]\n", "isoflux_set = np.array([[\n", " [4.25455147,\n", " 4.1881875,\n", " 4.2625,\n", " 4.45683769,\n", " 4.78746942,\n", " 5.31835938,\n", " 5.91875,\n", " 6.44275166,\n", " 6.92285156,\n", " 7.35441013,\n", " 7.73027344,\n", " 8.03046875,\n", " 8.19517497,\n", " 8.05673414,\n", " 7.75308013,\n", " 7.37358451,\n", " 6.94355469,\n", " 6.47773438,\n", " 5.99121094,\n", " 5.49433594,\n", " Rx,\n", " 4.79042969,\n", " 4.56269531,\n", " 4.36038615], \n", "[0.0,\n", " 1.13554688,\n", " 2.2565134,\n", " 3.16757813,\n", " 3.80507812,\n", " 4.0499021,\n", " 3.93089003,\n", " 3.665625,\n", " 3.31076604,\n", " 2.86875,\n", " 2.3199601,\n", " 1.62832118,\n", " 0.657421875,\n", " -0.35859375,\n", " -1.05585937,\n", " -1.59375,\n", " -2.03492727,\n", " -2.41356403,\n", " -2.75622016,\n", " -3.08699278,\n", " Zx,\n", " -2.40548044,\n", " -1.55982142,\n", " -0.67734375]\n", " ]])\n", " \n", "\n", "# instantiate the freegsnke constrain object\n", "from freegsnke.inverse import Inverse_optimizer\n", "constrain = Inverse_optimizer(null_points=null_points,\n", " isoflux_set=isoflux_set)\n", "\n", "\n", "# remove from the coils available for control the radial field coil \n", "# eq.tokamak['VS3'].control = False" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The inverse solve" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "GSStaticSolver.inverse_solve(eq=eq, \n", " profiles=profiles, \n", " constrain=constrain, \n", " target_relative_tolerance=1e-4,\n", " target_relative_psit_update=1e-3,\n", " # max_solving_iterations=150,\n", " # max_iter_per_update=10,\n", " # max_rel_update_size=0.075,\n", " # damping_factor=.99,\n", " # max_rel_psit=.05,\n", " verbose=True, # print output\n", " l2_reg=1e-14,\n", " )\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", "\n", "ax1.grid(zorder=0, alpha=0.75)\n", "ax1.set_aspect('equal')\n", "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", "constrain.plot(axis=ax1, show=False) # plots the contraints\n", "ax1.set_xlim(1.1, 12.4)\n", "ax1.set_ylim(-8.3, 8.3)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "eq.tokamak.getCurrents()\n", "\n", "# # save coil currents to file\n", "# import pickle\n", "# with open('simple_diverted_currents_PaxisIp.pk', 'wb') as f:\n", "# pickle.dump(obj=inverse_current_values, file=f)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "freegsnke4e", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.16" } }, "nbformat": 4, "nbformat_minor": 4 }